Climate-mediated shifts in temperature
fluctuations promote extinction risk

https://doi.org/10.1038/s41558-022-01490-7

Received: 2 March 2022

Kate Duffy

1,2,3

, Tarik C. Gouhier

4 & Auroop R. Ganguly1,5

Accepted: 31 August 2022

Published online: 20 October 2022

Check for updates

Climate-mediated changes in thermal stress can destabilize animal
populations and promote extinction risk. However, risk assessments
often focus on changes in mean temperatures and thus ignore the role of
temporal variability or structure. Using Earth System Model projections,
we show that significant regional differences in the statistical distribution
of temperature will emerge over time and give rise to shifts in the mean,
variability and persistence of thermal stress. Integrating these trends into
mathematical models that simulate the dynamical and cumulative effects
of thermal stress on the performance of 38 globally distributed ectotherm
species revealed complex regional changes in population stability over the
twenty-first century, with temperate species facing higher risk. Yet despite
their idiosyncratic effects on stability, projected temperatures universally
increased extinction risk. Overall, these results show that the effects of
climate change may be more extensive than previously predicted on the
basis of the statistical relationship between biological performance and
average temperature.

Biodiversity loss has been recognized as one of the top global risks
by the World Economic Forum because it could erode or eliminate
key ecosystem functions and services1. Climate change is expected to
surpass habitat loss as the leading threat to global biodiversity by the
middle of the twenty-first century2. Observed changes in the distri-
bution and phenology of species have already been linked to climate
fluctuations in numerous studies3. Although conservation actions
may ameliorate potential biodiversity loss, the success of these efforts
depends on our ability to predict the response of ecological systems
to environmental changes.

Most ecological impact studies so far have relied on statistical
models, such as bioclimate envelope approaches, to determine how cli-
mate change will impact ecological populations4–7. Bioclimate envelope
models are typically constructed by either mapping the geographical
distribution of species to co-located temperature records via regres-
sion techniques or by building species’ thermal profiles via empirical
assessments of their performance across a range of temperatures
(that is, thermal performance curves or TPCs)4,8. These relationships

between organisms and temperature are then used to predict the dis-
tribution of species under future thermal conditions projected under
various climate change scenarios.

Despite the power and popularity of TPCs, these statistical
approaches can yield inaccurate predictions because they typically
rely on mean annual conditions and thus ignore the influence of the
temporal structure of temperature fluctuations at finer scales. This
is problematic because the nonlinear relationship between tempera-
ture and most metrics of biological performance essentially guaran-
tees that the average organismal response will not be equivalent to
their response to the average condition9–12. Specifically, when an
organism is exposed to a sequence of temperatures x, its performance
at the average temperature f ( ̄x) will differ from the average of its
performance f (x) . Temporal variation in temperature will either
magnify (f (x) > f ( ̄x)) or dampen (f (x) < f ( ̄x)) the effects of its mean
on organismal performance depending on the curvature of f (that is,
whether f is accelerating or decelerating9). In many cases, changes in
temperature variability can be as or more important than changes in

1Sustainability and Data Sciences Laboratory, Department of Civil and Environmental Engineering, Northeastern University, Boston, MA, USA. 2NASA
Ames Research Center, Moffett Field, CA, USA. 3Bay Area Environmental Research Institute, Moffett Field, CA, USA. 4Department of Marine and
Environmental Sciences, Marine Science Center, Northeastern University, Nahant, MA, USA. 5Pacific Northwest National Laboratory, Richland, WA, USA.

e-mail: duffy.m.kate@gmail.com

Nature Climate Change | Volume 12 | November 2022 | 1037–1044

1037

nature climate change

Article

Page 2

the mean value13,14. In one study, climate-mediated changes in mean
temperature alone were found to broadly promote organismal per-
formance in ectotherms, but accounting for the temporal variability
of temperature dampened this effect and led to most species suffering
a performance loss15.

Although the temporal structure of temperature can theoretically
be incorporated into bioclimate envelope models by using finer tem-
poral scale data, accounting for its dynamical effects on organisms is
much more difficult because of the ‘static’ nature of these methods and
their general inability to account for the cumulative effects of previous
temperatures on organismal performance. However, theory has shown
that such carryover effects associated with the temporal structure
or autocorrelation of temperature can interact with the magnitude
of temperature variability to determine population persistence16.
Specifically, temporally autocorrelated variation tends to reduce
extinction risk by decreasing the likelihood of catastrophic conditions
under strong variation, whereas temporally autocorrelated variation
tends to promote extinction risk under weak variation by increasing
the likelihood that organisms will experience long stretches of poor
conditions16. Prolonged exposure to temperatures above the species’
critical thermal maximum is particularly destabilizing as it can reduce
population fitness below the replacement rate17. Analyses of historical
observations and projections from previous generation climate models
have found strong temporal trends in the variability and autocorrela-
tion of temperature18–21, suggesting the potential for a larger impact
on ecological populations in the future. Overall, these empirical and
theoretical results highlight the importance of quantifying changes
in the mean, variability and autocorrelation of temperature projected
under climate change to predict their joint influence on ecological
systems over the course of the twenty-first century. However, dispari-
ties in the scale of models in climate and ecology have hindered impact
studies that consider the complexity of both underlying systems22,23.
We briefly illustrate the potential for complex interactions
between climate-mediated changes in the mean, variability and auto-
correlation of temperature to influence organismal performance by
simulating the effects of synthetic temperature time series on the
population growth rate r according to a species’ TPC (Fig. 1, see Meth-
ods for modelling details). Predictably, performance under negligible
temperature variation can be inferred directly from the mean of each
species’ TPC (Fig. 1b,c). However, when temporal variation in tem-
perature is included in the model (that is, standard deviation; shaded
region), time-averaged performance can be considerably modified9,
even overturning the identity of ‘winning’ and ‘losing’ species based
solely on constant temperature conditions (Fig. 1d,e). Temperature
autocorrelation, which measures the temporal structure of tempera-
ture fluctuations (for example, the persistence of extremes), can also
play a pivotal role in determining whether a species’ performance
and stability will benefit or suffer under different thermal regimes
(Fig. 1f,g). To determine the impact of such changes over the course
of the twenty-first century, we analysed the latest generation of Earth
System Models from the Coupled Model Intercomparison Project
Phase 6 (CMIP6) to document spatiotemporal changes in three key
aspects of air temperature: statistical distribution, variance and tem-
poral autocorrelation. We then analysed the effects on population
stability and extinction risk using simple mathematical models to
examine the hypothesis that even under ideal conditions, popular
statistical methods can yield incorrect predictions about patterns of
organismal performance when dynamical and cumulative temperature
effects are ignored.

Regional trends in temperature distribution
We examined changes in the global and regional temperature distribu-
tions at each geographical location between 1850 and 2100 under the
high emissions scenario, SSP5-8.524 (Fig. 2a,b). Quantile regression was
used to measure temporal trends in the entire distribution of projected

temperatures (that is, across quantiles ranging from τ = 2.5% at the
low end to τ = 97.5% at the high end) in the Northern Hemisphere
Extra-tropics (NHEX, 30° N to 90° N), the Southern Hemisphere
Extra-tropics (SHEX, 90° S to 30° S), and the Tropics (TROP, 30° S to
30° N). When averaging trends across regions, we found asymmetrical
but uniformly positive trends across all quantiles, indicating that the
entire temperature distribution is shifting upwards but at rates that
vary systematically across the distribution. In NHEX, the lowest quan-
tile of the distribution (τ = 2.5%, 0.33 K per decade) is warming at twice
the rate of the uppermost quantile (τ = 97.5%, 0.16 K per decade). The
SHEX exhibits a similar pattern of disproportionate warming for the
low quantiles (τ = 2.5%, 0.15 K per decade; τ = 97.5%, 0.10 K per decade).
Conversely, in the tropics, the upper quantiles of temperature are warm-
ing faster (τ = 97.5%, 0.14 K per decade) than the lower quantiles
(τ = 2.5%, 0.10 K per decade). The magnitude of trends is greater in
NHEX than in SHEX or TROP. The more pronounced extra-tropical
decrease in the incidence of cold events may benefit cold-limited spe-
cies; however, quantile trends also indicate increased positive skewness
of the NHEX temperature distribution, which has been associated with
declines in long-term ecological performance15. Across all eight CMIP6
models that we analysed and in all three latitudinal regions, trends in
the tails of the distributions differ from the trends in the central ten-
dencies, thus highlighting the importance of moving beyond mean
temperature when predicting organismal performance.

Trends in the variability of temperature between 1850 and 2100
are predicted to exhibit similarly complex regional patterns (Fig. 2c).
Variance will generally increase across temperate and tropical land
areas below 45° N, with regional exceptions including Asia. The strong-
est increases in variance are in the northern mid latitudes, including
northern Africa, southern Europe, the Middle East and the western
United States. Variance is decreasing most rapidly in the high northern
latitudes, especially in Canada and Russia25. The concurrent decrease
in variability at high latitudes and its increase at other latitudes sug-
gests that temperature variation, similar to mean temperature, is
becoming more spatially homogeneous in a warming world. These
findings are generally consistent with studies of the previous genera-
tion of climate models, which suggested increasing temperature vari-
ability in tropical countries26 and decreasing variability in the northern
mid to high latitudes27. Trends at the regional level are congruent
with quantile trends (Fig. 2a), which indicate a widening temperature
distribution (increasing variance) in TROP, and a narrowing tempera-
ture distribution (decreasing variance) in NHEX and SHEX, as well as
large-scale changes in physical climate processes26–28. The effects
of these trends in temperature variation on ecological systems will
depend on the geographical location and physiological properties of
each species, with increasing variability either promoting or reducing
performance on the basis of its position relative to the inflection point
of an organism’s TPC9.

Frequency-resolved temperature changes
To better understand these spatiotemporal patterns, we used
time-frequency decomposition via the wavelet transform to resolve
changes in the variability of temperature at sub-annual to annual time-
scales (between 2 d and 2 yr) and multi-annual timescales (between 2 yr
and 30 yr; Extended Data Fig. 1). Wavelet transforms resolve a signal in
both the time and frequency domains to describe how each frequency
or period in the time series contributes to variation over time. We found
countervailing trends in scale-specific variability in the mid to high
northern latitudes. The magnitude of short-term variability is decreas-
ing, while the magnitude of long-term variability is increasing. Arctic
amplification, which is detectable in both observational data and cli-
mate simulations, has previously been suggested as the main driver of
decreasing sub-seasonal variability at these latitudes27. Meanwhile at
the mid latitudes, variation in both annual and multi-annual timescales
is increasing, consistent with increasing variance at all periodicities.

Nature Climate Change | Volume 12 | November 2022 | 1037–1044

1038

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 3

60° N

60° S

Cotesia sesamiae

Hyadaphis pseudobrassicae

60° E

120° E

180°

120° W

60° W

Mean of temperature

H. pseudobrassica
C. sesamiae

Variance of temperature

Temperature (°C)

)
g
o
l
(

r
e
w
o
P

Autocorrelation of temperature

β = –0.5
β = –2.0

Frequency (log)

30° N

0°

30° S

)
r
(

e
t
a
r

h
t
w
o
r
G

)

(

y
t
i
s
n
e
d
n
o
i
t
a
u
p
o
P

Time

Fig. 1 | Effects of temperature mean, variance and autocorrelation on
organismal performance. a, Source locations of the 38 species whose thermal
performance parameters were obtained from the Deutsch et al.5 dataset. Cotesia
sesamiae is a tropical parasitoid wasp and Hyadaphis pseudobrassicae is a
temperate-zone turnip aphid. b,c, Thermal performance curves and population
dynamics for C. sesamiae and H. pseudobrassicae under a mean temperature
(vertical line) with negligible variation. d,e, Larger temperature variation (s.d.,

shaded) alters mean response (dashed horizontal line) and may even overturn
predictions of relative performance based on constant temperature conditions.
f, The power spectrum of temperature with weak (ß = −0.5) and strong (ß = −2)
temporal autocorrelation. g, Population dynamics of Hyadaphis pseudobrassicae
under a greater degree of temporal autocorrelation exhibit longer-term
fluctuations. Multiple aspects of temperature, such as its mean and variance,
can interact to promote or decrease performance.

These scale-dependent changes in the temporal trends of temperature
fluctuations could have important ecological implications because the
effect of temperature fluctuations depends on the relationship between
their period and the generation time of organisms. Indeed, estimating
the biological effect of temperature fluctuations by ‘nonlinear averag-
ing’ of organismal performance under the relevant constant thermal

regimes is much more likely to yield accurate results when the period
of the temperature fluctuations is larger than the generation time of
an organism because such slow variation can more easily be ‘tracked’
by a population29.

We computed the spectral exponent of the temperature time
series at each geographical location to quantify spatiotemporal trends,

Nature Climate Change | Volume 12 | November 2022 | 1037–1044

1039

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 4

)
e
d
a
c
e
d

r
e
p
K

(

d
n
e
r
T

0.35

0.30

0.25

0.20

0.15

0.10

0.05

GLOBAL
NHEX
TROP
SHEX

30° N

0°

30° S

60° N

60° S

30° N

0°

30° S

60° N

60° S

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

120° W 60° W 0°

60° E 120° E

ΔK per decade

120° W 60° W 0°

60° E 120° E

ΔK 2 per decade

Temperature quantile

0.01

0.02

0.03

0.04

0.05

–0.5

–0.4

–0.3

–0.2

–0.1

0.1

Fig. 2 | Mean trends in the statistical distribution of daily air temperature
between 1850 and 2100. a,b, Trends in the percentile values of air temperature
(a, K per decade) and mean temperature at each geographic location (b, ΔK per
decade) indicate asymmetrically warming temperature distributions in the
Northern Hemisphere Extra-tropics (NHEX, 30° N to 90° N), the Tropics (TROP,
30° S to 30° N), the Southern Hemisphere Extra-tropics (SHEX, 90° S to 30° S),
and the full globe (GLOBAL, 90° S to 90° N). Shaded bounds denote the 90%

confidence interval based on eight CMIP6 models. c, Trends in the variance of
daily air temperature (ΔK2 per decade) exhibit similarly complex regional
patterns. The concurrent decrease in variability at high latitudes and increase at
other latitudes suggests that temperature variation is becoming more spatially
homogeneous in a warming world. Hashed contours indicate statistically
significant inter-model agreement on the sign of the trend at the α = 0.05
significance level.

with more negative exponents indicating greater temporal autocorrela-
tion over a range of lags from 2 d to 10 yr (Fig. 3a). We found increasing
temporal autocorrelation (decreasing spectral exponent) at a majority
of sea locations (60%) and land locations (80%), excluding Antarctica
where autocorrelation is decreasing. Autocorrelation is increasing
most rapidly in equatorial land areas including the Amazon and the
Southeast Asian islands, with high inter-model agreement on the sign
of the trend. Notable exceptions to the increasing trend in autocorrela-
tion include Greenland, Western Africa, Western Europe and parts of
Central Asia. Generally, agreement between models is higher at mid lati-
tudes than in the polar zones or the tropics, where climate model bias
and spread have historically persisted30. Regional analysis indicates
statistically significant increasing trends in temporal autocorrelation
in NHEX (−1.12 × 10−3 per decade, P = 0.010), TROP (−1.14 × 10−3 per dec-
ade, P = 0.001) and globally (−0.54 × 10−3 per decade1, P = 0.005), and a
statistically significant decreasing trend in temporal autocorrelation
in SHEX (0.53 × 10−3 per decade, P = 0.009; Supplementary Table 1).
The direction and significance of these trends are consistent across
land and sea environments, although the spectral exponent is more
negative for sea than land, probably due to the buffering effects of the
ocean (Fig. 3b–e). In NHEX and TROP, autocorrelation is increasing at
a greater rate in land locations than in sea locations, while in SHEX it
is decreasing at similar rates between land and sea (Supplementary
Table 2). A greater degree of temporal autocorrelation is associated
with more gradual changes of state and, even without any changes in
variance, results in longer durations spent under extreme conditions.
A greater clustering of similar temperatures has been suggested to
increase exposure to heat waves and cold snaps while decreasing the
incidence of protective temporal refugia20.

Regional differences in warming patterns
In the northern latitudes, variance and autocorrelation exhibit oppo-
site temporal trends. The decreasing variance may be attributed to
a decrease in high-frequency variability and more rapid warming of
the lower than the upper quantiles of the temperature distribution.
Studies of reanalysis data and observations have also implicated
decreasing cold-season sub-seasonal variability and rapidly warming
cold days in decreasing temperature variability in mid to high north-
ern latitudes20,24,29. Meanwhile, temporal autocorrelation in NHEX is
increasing—a finding that has also been detected in the previous genera-
tion of climate models20, weather station observations31 and monthly
reanalysis data19. As a result, variation at 2 d to 10 yr periodicities is
decreasing while temperature fluctuations are becoming more persis-
tent, suggesting the increased probability of a series of homogeneous

conditions. In contrast to the mid to high northern latitudes, variance
and temporal autocorrelation show similar trends at most latitudes,
that is, both variance and autocorrelation are increasing.

Implications for global ectotherm populations
To better understand the independent and joint effects of these pro-
jected trends in the mean, variance and autocorrelation of temperature
on ecological systems, we used empirical thermal performance infor-
mation from invertebrate ectotherms compiled by Deutsch et al.5. We
extracted temperature time series from the eight CMIP6 climate mod-
els at geographical point locations corresponding to the source sites
of the 38 species (Fig. 4a). A dynamical population simulation using
species-specific temperature-dependent growth rates yielded time
series of population abundance for the historical period (1950–2000)
and the latter half of the twenty-first century (2050–2100). We used a
dynamical logistic growth model whose carrying capacity K = rt/α is
determined by the temperature-dependent growth rate rt and the
self-regulation parameter α. Importantly, the model captures the
effects of temperatures above the critical thermal maximum and extinc-
tion propensity under autocorrelated variation by allowing growth
rates to become negative (see Methods for details). Using the eight
climate simulations as replicates, we compared the historical and future
periods to detect statistically significant temperature-driven changes
in population abundance, stability (mean/standard deviation of abun-
dance) and extinction probability (proportion of simulations where a
species did not have a strictly positive final abundance).

Under the high emissions scenario (SSP5-8.5), population abun-
dance increased for the plurality of species (16 of 38) because the mean
temperature grew closer to their thermal optimum and thus boosted
equilibrium abundance, but it decreased for 9 species (Supplementary
Table 3). Population abundance increased significantly for 3 of 5 TROP
species and for the majority (5 of 8) of SHEX species. In NHEX, outcomes
were mixed, with approximately equal proportions of species experi-
encing an increase in abundance, a decrease in abundance, and no
significant change. NHEX population abundance followed latitudinal
patterns, generally decreasing between 30° N and 45° N, and increasing
above of 45° N. Under the high emissions scenario, population stability
increased for the 12 out of 38 and decreased for 9 species (Fig. 4b).
Population stability increased or underwent no significant change for
TROP species, while in the mid latitudes (NHEX and SHEX), changes in
stability were mixed. Additional analyses showed that the trends in
stability were mainly due to the emergence of two distinct dynamical
regimes under climate change, with species either moving to a
low-mean/low-variance mode or a high-mean/high-variance mode,

Nature Climate Change | Volume 12 | November 2022 | 1037–1044

1040

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 5

30° N

0°

30° S

60° N

60° S

120° W 60° W

0°

60° E 120° E

Change in spectral exponent per decade

–1.0

–0.5

0.5

Fig. 3 | Increasing temporal autocorrelation in daily air temperature
between 1850 and 2100. a, Spatiotemporal trends in temporal autocorrelation
suggest changes in the chronological sequence of temperature conditions, with
increasing temporal autocorrelation (decreasing spectral exponent) at 80.04%
of global land locations, excluding Antarctica. Hashed contours indicate
statistically significant inter-model agreement on the sign of the trend at the
α = 0.05 significance level. b–e, Regional analysis indicates statistically

t
e
n
o
p
x
e

a
r
t
c
e
p
S

t
e
n
o
p
x
e

a
r
t
c
e
p
S

–1.24

–1.26

–1.28

–1.30

–1.32

–1.34

–1.24

–1.26

–1.28

–1.30

–1.32

–1.34

GLOBAL

Land
Sea

1850

1900

1950

2000

2050

2100

NHEX

Land
Sea

t
e
n
o
p
x
e

a
r
t
c
e
p
S

t
e
n
o
p
x
e

a
r
t
c
e
p
S

–1.36

–1.38

–1.40

–1.42

–1.44

–1.46

–1.12

–1.14

–1.16

–1.18

–1.20

–1.22

–1.24

TROP

Land
Sea

1850

1900

1950

2000

2050

2100

SHEX

Land
Sea

1850

1900

1950

2000

2050

2100

1850

1900

1950

2000

2050

2100

significant increasing trends in temporal autocorrelation in NHEX and TROP,
and a statistically significant decreasing trend in temporal autocorrelation in
SHEX. While sea environments generally exhibit a greater degree of temporal
autocorrelation than land, in NHEX autocorrelation is increasing at a greater
rate on land locations as to overturn this relationship by the end of the
twenty-first century.

particularly in the extra-tropics (Extended Data Figs. 2 and 3). These
results were robust to orders of magnitude changes in the growth
rate rt and self-regulation parameter α (Extended Data Figs. 4 and 5).
Many SHEX and NHEX species suffered performance losses (nega-
tive growth rates) during summers in their respective hemispheres, as
they are generally less tolerant of hot temperatures than tropical spe-
cies. For some temperate species, longer growing seasons and warmer
winter temperatures offset the negative effect of the warmest part of
the year, while others suffered an overall performance loss32. This is
consistent with the suggestion that increases in summer heat stress
would reduce overall fitness and increase fitness variation for many
mid-latitude species. Our results suggest that temperate species may
be at greater risk than tropical species as a result of warm days, even
when annual mean temperature remains below the thermal optimum.
The results contrast with those of previous studies, which suggested
on the basis of hourly temperature records and monthly temperature
anomalies that warming in the tropics would be more deleterious than
warming in the mid latitudes5,33. This discrepancy may be due to the fact
that growth rates were allowed to become negative when temperatures
exceeded the critical thermal maximum in our simulations but assumed
to converge to zero (that is, were not allowed to be negative) in previous
studies4. Our results are more consistent with studies that predict a
greater risk of performance loss for temperate species when account-
ing for negative performance values in response to climate-mediated
changes in the mean and the variance of temperature15.

To tease apart the dynamical effects of climate change on popu-
lation stability from its effects on mean performance as inferred by
measuring average growth rate using each species’ TPC, we replicated
previous efforts by comparing changes in the average growth rate
under historical and future climatic conditions with vs without negative
growth rates (Extended Data Fig. 6). Our results show that although
allowing negative growth rates predictably leads to greater reduc-
tions in performance overall, the regional patterns in performance are
similar to the trends in population stability observed in the dynamical
simulations, with tropical species generally enjoying performance
gains and temperate species—particularly in NHEX—suffering perfor-
mance losses (Extended Data Fig. 6).

Our simulations indicated mean warming as the dominant driver
of ecological impacts. Changes in temporal autocorrelation alone
(mean temperature and variance held at historical levels) had no

significant effects on population abundance and a significant desta-
bilizing effect on just 1 NHEX species. Changes in temporal autocor-
relation and variance (mean temperature held at historical levels)
led to an increase in population abundance in 2 NHEX species and a
decrease in population stability in 2 NHEX species. These results sug-
gest that NHEX species are more susceptible to changes in temperature
variability than TROP or SHEX species. Finally, changes in mean and
temporal autocorrelation (variance held at historical levels) led to
increased population abundance in 18 global species and increased
stability in 14 global species, versus 16 and 12 under the high emissions
scenario-projected changes in all three aspects of temperature. Thus,
projected changes in temperature variability have a weak moderating
effect on the positive effects of mean warming on population abun-
dance and stability.

To determine how these complex changes in population abun-
dance and stability translate to persistence, we quantified extinction
risk as the proportion of the 8 CMIP6 models for which population
abundance declined below an arbitrarily small threshold of 1 × 10−9 at
any point during the 50 yr simulation (Fig. 4c). In our simulations under
the high emissions scenario, extinction risk increased significantly
under future climate conditions relative to historical baselines for
18 species, increased (but not significantly) for 6 species, decreased
for 1 species, and did not change for 13 species. We found statisti-
cally significant increases in extinction risk globally (Mann–Whitney
U = 423, n1 = n2 = 38, P = 8 × 10−4) and in NHEX (Mann–Whitney U = 166,
n1 = n2 = 25, P = 3 × 10−3). These findings suggest that temperature
changes promote extinction risk, despite having a largely positive or
neutral effect on population abundance and idiosyncratic impacts on
stability. Hence, although variability among climate models produces
a wide range of changes in stability across species and geographical
locations, uncertainty at the climate level yields consistent biological
impacts in the form of systematically higher extinction risks (Extended
Data Fig. 7).

Conclusion
By forcing simple strategic and dynamical models of population growth
with fine temporal scale temperature projections from the latest gen-
eration of Earth System Models, we demonstrated increased extinc-
tion risk under climate change across globally distributed ectotherm
populations. Unfortunately, using more complex tactical dynamical

Nature Climate Change | Volume 12 | November 2022 | 1037–1044

1041

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 6

y
t
i
l
i

b
a
t
s

n
o
i
t
a
u
p
o
p
n

e
g
n
a
h
c

t
n
e
c
r
e
P

60° N

30° N

0°

30° S

60° S

NHEX
TROP
SHEX

300

200

100

–100

300

200

100

–100

300

200

100

–100

300

200

100

–100

23
22

35
34

32
31 30
28
29

18
16

9
6
5
3

60° E

120° E

180°

120° W

60° W

Mean, variance and autocorrelation

Mean and autocorrelation

Variance and autocorrelation

Autocorrelation

y
t
i
l
i

b
a
b
o
r
p
n
o
i
t
c
n
i
t
x
e

e
g
n
a
h
C

100

–50

100

–50

100

–50

100

–50

Mean, variance and autocorrelation

NHEX
TROP
SHEX

Mean and autocorrelation

Variance and autocorrelation

Autocorrelation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

Species

Fig. 4 | Temperature has idiosyncratic effects on stability but increases
extinction risk globally. a, Source locations of the terrestrial ectothermic
invertebrate species, numbered 1 (southernmost latitude) to 38 (northernmost
latitude). Species are colour-coded according to latitudinal region (orange,
SHEX; green, TROP; red, NHEX). b, Percent changes in population stability
(mean ÷ s.d.) between a historical reference period (1950–2000) and a future
period (2050–2100) under multiple aspects of temperature change indicate

greater risk to temperate than to tropical species. Under a high emissions
scenario, stability shows a statistically significant increase for 12 of 38 species
and a statistically significant decrease for 9 species. Points in the violin plots
represent the 8 climate model outputs. c, Extinction probability shows a
quasi-universal increase globally between the historical period (1950–2000)
and a future period (2050–2100) under high emissions scenario changes in
temperature.

models would require extensive species-, age- and life-stage-specific
information about the effects of temperature fluctuations on popu-
lation growth rates that is simply not available at the relevant scales.
Tactical models would also need to consider thermoregulation34,
the effects of microclimates35, acclimatization or adaptation36, par-
titioning of activity periods37 and synecological processes such as
predator-prey interactions that could affect ectotherm population

dynamics. Additionally, due to their 1° spatial resolution, the climate
projections used in this study are much coarser than the microclimates
experienced by individual organisms and may thus lead to underes-
timates of organismal performance due to the presence of thermal
refugia in the real world23,34. Hence, our results should be viewed as a
qualitative baseline prediction of how the spatiotemporal distribution
of extinction risk is likely to shift due to climate change rather than a

Nature Climate Change | Volume 12 | November 2022 | 1037–1044

1042

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 7

quantitative forecast of when each species is likely to be extirpated
from each geographical location.

Despite the limitations of TPCs in accounting for temporal car-
ryover and dynamical effects, the lack of obvious alternatives calls for
strategies to make these approaches more robust to real-world condi-
tions38, such as by integrating more realistic, fine-scaled temperature
variation into our predictive models than previous studies. Although
bioclimate envelope approaches have been criticized for not account-
ing for important ecological factors, such as species interactions and
dispersal, when attempting to predict the ecological effects of climate
change39–42, we have shown that even under ideal conditions when the
influence of such factors can be assumed to be negligible, statistical
frameworks that ignore the dynamical consequences of temperature
variation are likely to yield inaccurate forecasts of the impact of climate
change on organisms. Our results show that accounting for shifts in the
entire statistical distribution of temperature over time via dynamical
models can better capture the cumulative effects of climate-mediated
changes in thermal stress on extinction risk.

By bringing together climate data and a minimal dynamical model
from ecology, we demonstrated a strong and systematic amplifica-
tion of extinction risk in ectotherms due to projected changes in
fine-grained temperature variability. Furthermore, our finding of
greater risk to sub-tropical than tropical species highlights the impor-
tance of accounting for the dynamical effects of projected changes in
the mean as well as the variance of temperature over the course of the
twenty-first century to accurately predict the response of ecological
systems around the globe.

Online content
Any methods, additional references, Nature Research reporting sum-
maries, source data, extended data, supplementary information,
acknowledgements, peer review information; details of author contri-
butions and competing interests; and statements of data and code avail-
ability are available at https://doi.org/10.1038/s41558-022-01490-7.

References
1.

The Global Risks Report 2021 Insight Report 16th edn (World
Economic Forum, 2021).

2. Bellard, C., Bertelsmeier, C., Leadley, P., Thuiller, W. & Courchamp,
F. Impacts of climate change on the future of biodiversity. Ecol.
Lett. 15, 365–377 (2012).

12. Denny, M. Performance in a variable world: using Jensen’s

inequality to scale up from individuals to populations. Conserv.
Physiol. 7, coz053 (2019).

13. García-Carreras, B. & Reuman, D. C. Are changes in the mean
or variability of climate signals more important for long-term
stochastic growth rate? PLoS ONE 8, e63974 (2013).

14. Benedetti-Cecchi, L., Bertocci, I., Vaselli, S. & Maggi, E. Temporal

variance reverses the impact of high mean intensity of stress in
climate change experiments. Ecology 87, 2489–2499 (2006).
15. Vasseur, D. A. et al. Increased temperature variation poses a

greater risk to species than climate warming. Proc. R. Soc. B 281,
20132612 (2014).

16. Schwager, M., Johst, K. & Jeltsch, F. Does red noise increase or

decrease extinction risk? Single extreme events versus series of
unfavorable conditions. Am. Nat. 167, 879–888 (2006).
17. Kingsolver, J. G., Diamond, S. E. & Buckley, L. B. Heat stress

and the fitness consequences of climate change for terrestrial
ectotherms. Funct. Ecol. 27, 1415–1423 (2013).

18. Meehl, G. A. More Intense, more frequent, and longer

lasting heat waves in the 21st century. Science 305,
994–997 (2004).

19. Lenton, T. M., Dakos, V., Bathiany, S. & Scheffer, M. Observed

trends in the magnitude and persistence of monthly temperature
variability. Sci. Rep. 7, 5940 (2017).

20. Di Cecco, G. J. & Gouhier, T. C. Increased spatial and temporal

autocorrelation of temperature under climate change. Sci. Rep. 8,
14850 (2018).

21. Li, J. & Thompson, D. W. J. Widespread changes in surface

temperature persistence under climate change. Nature 599,
425–430 (2021).

22. Lembrechts, J. J. et al. Comparing temperature data sources for

use in species distribution models: from in-situ logging to remote
sensing. Glob. Ecol. Biogeogr. 28, 1578–1596 (2019).

23. Potter, K. A., Arthur Woods, H. & Pincebourde, S. Microclimatic
challenges in global change biology. Glob. Change Biol. 19,
2932–2939 (2013).

24. O’Neill, B. C. et al. The Scenario Model Intercomparison
Project (ScenarioMIP) for CMIP6. Geosci. Model Dev. 9,
3461–3482 (2016).

25. Hansen, J., Sato, M. & Ruedy, R. Perception of climate change.

Proc. Natl Acad. Sci. USA 109, E2415–E2423 (2012).

3. Parmesan, C. Ecological and evolutionary responses to

26. Bathiany, S., Dakos, V., Scheffer, M. & Lenton, T. M. Climate models

recent climate change. Annu. Rev. Ecol. Evol. Syst. 37,
637–669 (2006).

4. Pearson, R. G. & Dawson, T. P. Predicting the impacts of climate
change on the distribution of species: are bioclimate envelope
models useful? Glob. Ecol. Biogeogr. 12, 361–371 (2003).
5. Deutsch, C. A. et al. Impacts of climate warming on terrestrial
ectotherms across latitude. Proc. Natl Acad. Sci. USA 105,
6668–6672 (2008).

predict increasing temperature variability in poor countries. Sci.
Adv. 4, eaar5809 (2018).

27. Screen, J. A. Arctic amplification decreases temperature variance

in northern mid- to high-latitudes. Nat. Clim. Change 4, 577–582
(2014).

28. Stouffer, R. J. & Wetherald, R. T. Changes of variability in response

to increasing greenhouse gases. Part I: temperature. J. Clim. 20,
5455–5467 (2007).

6. Cheung, W. W. L. et al. Projecting global marine biodiversity

29. Gouhier, T. C. & Pillai, P. Commentary: nonlinear averaging

impacts under climate change scenarios. Fish Fish. 10, 235–251
(2009).
Thuiller, W. et al. Consequences of climate change on the tree of
life in Europe. Nature 470, 531–534 (2011).

8. Angilletta, M. J. Thermal Adaptation: A Theoretical and Empirical
Synthesis. Thermal Adaptation (Oxford Univ. Press, 2009).
9. Ruel, J. J. & Ayres, M. P. Jensen’s inequality predicts effects of
environmental variation. Trends Ecol. Evol. 14, 361–366 (1999).

of thermal experience predicts population growth rates in a
thermally variable environment. Front. Ecol. Evol. 7, 236 (2019).
30. Tian, B. & Dong, X. The double-ITCZ bias in CMIP3, CMIP5, and
CMIP6 models based on annual mean precipitation. Geophys.
Res. Lett. 47, e2020GL087232 (2020).

31. Dillon, M. E. et al. Life in the frequency domain: the biological

impacts of changes in climate variability at multiple time scales.
Integr. Comp. Biol. 56, 14–30 (2016).

10. Lawson, C. R., Vindenes, Y., Bailey, L. & Pol, Mvande

32. Adamo, S. A., Baker, J. L., Lovett, M. M. E. & Wilson, G.

Environmental variation and population responses to global
change. Ecol. Lett. 18, 724–736 (2015).

11. Denny, M. The fallacy of the average: on the ubiquity, utility

Climate change and temperate zone insects: the tyranny of
thermodynamics meets the world of limited resources. Environ.
Entomol. 41, 1644–1652 (2012).

and continuing novelty of Jensen’s inequality. J. Exp. Biol. 220,
139–146 (2017).

33. Dillon, M. E., Wang, G. & Huey, R. B. Global metabolic impacts of

recent climate warming. Nature 467, 704–706 (2010).

Nature Climate Change | Volume 12 | November 2022 | 1037–1044

1043

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 8

34. Sunday, J. M. et al. Thermal-safety margins and the necessity of
thermoregulatory behavior across latitude and elevation. Proc.
Natl Acad. Sci. USA 111, 5610–5615 (2014).

35. Pincebourde, S. & Casas, J. Narrow safety margin in the

phyllosphere during thermal extremes. Proc. Natl Acad. Sci. USA
116, 5588–5596 (2019).

36. Somero, G. N. The physiology of climate change: how potentials
for acclimatization and genetic adaptation will determine
‘winners’ and ‘losers’. J. Exp. Biol. 213, 912–920 (2010).

37. Johansson, F., Orizaola, G. & Nilsson-Örtman, V. Temperate insects
with narrow seasonal activity periods can be as vulnerable to
climate change as tropical insect species. Sci. Rep. 10, 8822 (2020).

38. Sinclair, B. J. et al. Can we predict ectotherm responses to

climate change using thermal performance curves and body
temperatures? Ecol. Lett. 19, 1372–1385 (2016).

39. Davis, A. J., Jenkinson, L. S., Lawton, J. H., Shorrocks, B. & Wood,
S. Making mistakes when predicting shifts in species range in
response to global warming. Nature 391, 783–786 (1998).

40. Suttle, K. B., Thomsen, M. A. & Power, M. E. Species interactions
reverse grassland responses to changing climate. Science 315,
640–642 (2007).

42. Harley, C. D. G. Climate change, keystone predation, and

biodiversity loss. Science 334, 1124–1127 (2011).

Publisher’s note Springer Nature remains neutral with
regard to jurisdictional claims in published maps and
institutional affiliations.

Open Access This article is licensed under a Creative Commons
Attribution 4.0 International License, which permits use, sharing,
adaptation, distribution and reproduction in any medium or format,
as long as you give appropriate credit to the original author(s) and the
source, provide a link to the Creative Commons license, and indicate
if changes were made. The images or other third party material in this
article are included in the article’s Creative Commons license, unless
indicated otherwise in a credit line to the material. If material is not
included in the article’s Creative Commons license and your intended
use is not permitted by statutory regulation or exceeds the permitted
use, you will need to obtain permission directly from the copyright
holder. To view a copy of this license, visit http://creativecommons.
org/licenses/by/4.0/.

41. Gouhier, T. C., Guichard, F. & Menge, B. A. Ecological processes
can synchronize marine population dynamics over continental
scales. Proc. Natl Acad. Sci. USA 107, 8281–8286 (2010).

© This is a U.S. Government work and not under copyright protection
in the US; foreign copyright protection may apply 2022, corrected
publication 2023

Nature Climate Change | Volume 12 | November 2022 | 1037–1044

1044

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 9

Methods
CMIP6 simulations
We obtained CMIP6 climate simulations for the historical forcing
period (1850–2014) and future emissions scenario SSP5-8.5 (2015–
2100) via the CMIP6 data portal (https://esgf-node.llnl.gov/search/
cmip6/). Eight models from CMIP6 (AWI-CM-1-1-MR, BCC-CSM2-MR,
CESM2, EC-Earth3, INM-CM5-0, MPI-ESM1-2-HR, MRI-ESM2-0 and
NorESM2-MM) were selected on the basis of availability of daily air
temperature at surface (‘tas’) at a 100 km nominal resolution at the
time of download. While ‘tas’ at sub-daily frequencies is available
for some models, daily data was selected to maximize the ensemble
size. We resampled all datasets to a common 1° by 1° grid spanning
−90° to 90° latitude and 0° to 360° longitude, and to a standard cal-
endar without leap years. Spatial regions were defined on the basis
of latitude as Northern Hemisphere Extra-tropics (90° S to 30° S),
Tropics (30° S to 30° N) and Southern Hemisphere Extra-tropics
(30° N to 90° N).

Statistical analyses of climate data
Quantile regression. Trends in the percentile values of global and
regional temperature distributions were computed via quantile
regression. Quantile regression can comprehensively model hetero-
geneous conditional distributions, where the relationship between
the quantiles of the dependent variable and the independent variable
is different from the relationship between the means of the depend-
ent variable and the independent variable. We applied quantile
regression to analyse trends with respect to time at various percentile
values (P2.5, P10, P20, P30, P40, P50, P60, P70, P80, P90, P97.5). Analyses were
performed using the R package quantreg, with significance level
α = 0.1 and the default Barrodale and Roberts method to return
confidence intervals for the estimated parameters. To obtain the
ensemble mean trends, we calculated the mean slope, upper bound
and lower bound across the eight climate models at each geographi-
cal location, then computed spatial averages for the full globe and
three latitudinal regions.

Variance. Trends in the magnitude of temporal variation of air tem-
perature were examined at each geographical location using a mov-
ing window approach. First, temperature was detrended by fitting a
piecewise linear regression against time with Python package pwlf
at each geographical location and extracting the residuals. Then, the
temperature time series were divided into 10 yr windows starting in
years 1855 through 2085 so as not to combine historical and future
simulations (pre- and post-01-01-2015), and the variance of daily air
temperature was calculated for each window. Windows were selected
with no overlap to avoid statistical issues due to non-independence of
estimates taken from partially overlapping time windows20.

Scale-specific variability. Scale-specific variability was quantified
using time-frequency decomposition. Specifically, at each geographi-
cal location, wavelet analysis was conducted on multi-model mean
temperature using the R package biwavelet43. Wavelet analysis resolves
both the time and frequency domains of a signal (here a time series) via
the wavelet transform. This is achieved via the convolution of a mother
wavelet function and a time series across a set of windows τ and scales s.
We chose to use the Morlet wavelet, which represents a sine wave
modulated by a Gaussian function:44

ψ0 (t) = π−1/4eiωt

0 e−r2/2

where i is the imaginary unit, t represents non-dimensional time, and
ω0 = 6 is the non-dimensional frequency3. The continuous wavelet
transform of a discrete time series x (t) with equal spacing δt and length
T is defined as the convolution of x (t) with a normalized Morlet
wavelet:44,45

Nature Climate Change

Wx (s, τ) =

√√√
√

δt
s

T−1
∑
t=0

x (t) ψ0 ∗ (

(t − τ) δt
s

)

where * indicates the complex conjugate. By varying the wavelet scale
s (that is, dilating and contracting the wavelet) and translating along
localized time position τ, one can calculate the wavelet coefficients
Wx (s, τ) across the different scales s and positions τ. These wavelet coef-
ficients can be used to compute the bias-corrected local wavelet power,
which describes how the contribution of each frequency or period in
the time series varies over time:44,46,47

x (s, τ) = 2s |Wx (s, τ)|

2where 2s is the bias correction factor46. The

scale s of the Morlet wavelet is related to the Fourier frequency f:47,48

4πs
ω0+√2+ω2

When ω0 = 6, the scale s is approximately equal to the reciprocal
of the Fourier frequency f, so period p ≈ s. The local wavelet power
spectrum can then be visualized via heat maps and contour plots45,47.
From the resulting local wavelet power spectrum heat map with time
on the x axis, period (scale) on the y axis and power on the z axis,
scale-averaged wavelet power was computed at annual (between 3 d
and 2 yr) and multi-annual (between 2 yr and 30 yr) periodicities. This
was achieved by taking the weighted sum of the local wavelet power
across all scales for each time location τ:44,47

W−2
x

(τ) =

δjδt
Cδ

J
∑
j=0

2
||Wx (sj, τ)||
sj

where Cδ = 0.776 for the Morlet wavelet, δJ represents the spacing
between successive scales and δt represents the spacing between suc-
cessive time locations44. Scale-averaged power was then regressed
against time using Generalized Least Squares (GLS) regression for the
period 1850–2100 at each geographic location. To determine the
robustness of results to the choice of period for scale averaging, we
also performed analysis of trends separately at interannual (between
2 yr and 7 yr) and multi-annual (between 7 yr and 30 yr) scales and found
qualitatively similar results.

Temporal autocorrelation. The temporal autocorrelation of air
temperature was quantified by calculating the spectral exponent at
each geographical location20. As described above, temperature was
detrended by fitting a piecewise linear regression at each geographical
location and extracting the residuals. The detrended temperature was
divided into 10 yr windows starting in years 1855 through 2085. Fourier
transforms of each time series were computed via fast Fourier trans-
form using the Python package NumPy. Periodograms were prepared
with frequency on the x axis and power spectral density on the y axis.
The spectral exponent, β, was calculated as the slope of the regression
line relating log-transformed power to log-transformed frequency.
β expresses the relative contributions of frequencies to the power
spectrum. In the case of equal contribution from all frequencies, β = 0.
Greater contribution from low frequencies than from high frequencies
results in a more negative value of β and indicates greater temporal
autocorrelation in the time domain.

Analysis of decadal trends. For each climate model, GLS regression
was used to detect statistically significant trends (P < 0.05) in variance
and temporal autocorrelation with respect to time in the presence of
potentially autocorrelated residuals. To measure inter-model agree-
ment, we calculated the multi-model mean trend as the mean of trends
calculated for each of the 8 models at each geographic location, then
assessed the proportion of models that agreed with the sign of the
multi-model mean trend. Inter-model agreement was considered as
statistically significant at the α = 0.1 level on the basis of a binomial
test. ANCOVA was used to quantify the relationship between temporal

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 10

autocorrelation and time while accounting for potential differences
between land and sea environments. Statistically significant main
effects and interactions were reported for P < 0.05.

dynamics. Specifically, we used the r − α logistic growth model to
simulate temperature-dependent growth and negative
density-dependence:

Modelling temperature impacts on ecology
Thermal tolerance data. We obtained experimentally derived ther-
mal tolerance parameters for a set of terrestrial ectotherms (n = 38)
published by Deutsch et al.5 and used them to predict physiological
response to CMIP6-simulated temperature. Deutsch et al.5 gathered
data from 31 thermal performance studies published between 1974
and 2003 based on a collection of insects from 35 different locations.
For each species, experimental intrinsic growth rates at multiple tem-
peratures were used to fit a TPC, yielding least-squares estimates of
key parameters such as critical thermal maximum (CTmax), optimum
temperature (Topt), and sigma (σ). We used a numerical scheme to
reconstruct the curves whereby the rise in performance up to Topt was
modelled as Gaussian and the decline beyond Topt was quadratic5,49:

P (T) =

⎧
⎪

⎨
⎪
⎩

exp [− (

T−Topt

2σ

2
)

] for T ≤ Topt

1 − (

T−Topt

Topt−CTmax

)

for T > Topt

(1)

This allowed negative growth rates to arise at high temperatures,
but growth rates were bound at zero at low temperatures. Negative
performance values indicate that mortality surpasses reproduction
rates. Because P(T) is capped at 1 under this numerical scheme, P(T)
represents the relative fitness of each species based on its normalized
maximum growth rate. However, scaling this relative or normalized
maximum growth rate by two orders of magnitude (that is, by a factor
of 0.1 or 10.0) had limited quantitative and no qualitative impact on
our results (Extended Data Fig. 4). Overall, increasing the growth rate
scaling factor had no impact on population stability but promoted
extinction risk.

Isolation of temperature aspects. To isolate projected changes in
mean temperature and variability, we transformed the future (2050–
2100) time series using z-score normalization. Using this approach, we
modified projected time series to match the historical (1950–2000)
mean and/or standard deviation. Working in 10 yr moving windows
between 2050 and 2100, each series xi with mean m1 and standard
deviation s1 was transformed to series yi with mean m2 and standard
deviation s2:

yi = m2 + (xi − m1)

s2
s1

(2)

According to the scenario, m2 and s2 were alternatively defined
as (1) high emissions scenario mean and standard deviation (‘Mean,
variance and autocorrelation’), (2) high emissions scenario mean and
historical standard deviation (‘Mean and autocorrelation’), (3) histori-
cal mean and high emissions scenario standard deviation (‘Variance
and autocorrelation’) and (4) historical mean and standard deviation
(‘Autocorrelation’). High emissions scenario statistics refer to the
properties of future series xi and confer no change to that aspect of
the time series.

Population dynamical modelling. To model the effects of tempera-
ture change on the stability and extinction probability of global ecto-
therm populations, temperature dependence was integrated in the
growth rate term of a population dynamical model50. While more com-
plex synecological models can capture a range of community-level
effects including competition and predation, we chose to model first
order autecological dynamics to produce foundational insights about
the role of temperature fluctuations on single-species population

Nature Climate Change

dN
dt

= N (rt − αN)

(3)

with population size N, time t, temperature-dependent growth rate rt,
and self-regulation in the form of intraspecific competition α. This r − α
logistic model is easily interconvertible with the classical r − K formula-
tion (r/α = K), but has the advantage of handling negative values of r
without issues51. This approach is sensitive to the effects of tempera-
tures at and above the critical thermal maximum, which can yield
negative growth rates that are important for determining population
dynamics as well as long-term fitness.

We extracted time series of daily temperature at the source loca-
tions for each species from the ensemble of eight climate simulations.
Daily intrinsic growth rates were computed from temperature using
equation (1), incorporated into the r − α logistic growth model depicted
in equation (3), and the model was then numerically solved using the
explicit Runge-Kutta method of order 5(4) implemented in the Python
SciPy package to obtain daily population densities. Rather than delin-
eating active periods, which may shift under climate change, we con-
sidered the full year to account for potential changes in fitness due to
shifts in activity.

The sensitivity of the results to strong (α = 1) and weak (α = 0.1)
self-regulation was examined and found to be extremely limited
(Extended Data Fig. 5). We also assessed the sensitivity of our results to
absolute rather than relative or normalized growth rates by scaling rt by
a factor of 0.1 or 10 in our simulations. Scaling rt by two orders of mag-
nitude in this manner had very little quantitative and no qualitive impact
on our results. This suggests that the effects of temperature fluctuations
on changes in the spatiotemporal distribution of population abundance,
stability and extinction were not contingent upon the use of relative
fitness (that is, normalized growth rate) versus absolute fitness (that is,
growth rate scaled by a factor of 0.1 or 10). These sensitivity analyses also
showed that our results are robust to temperature-mediated changes in
the maximum instantaneous growth rate52,53.

Analysis of population changes. To quantify temperature-driven
changes in ecological stability and extinction probability, we compared
population sizes and dynamics between a historical period (1950–
2000) and a future period (2050–2100). Here we defined latitudinal
regions according to traditional delineations in ecology: Northern
Hemisphere Extra-tropics, 60° S to 23° S; Tropics, 23° S to 23° N; and
Southern Hemisphere Extra-tropics, 23° N to 60° N.

Population abundance was computed as the mean population
size (N) for a time period. Population stability was computed as the
inverse of the coefficient of variation, or mean population divided
by population standard deviation. Percent changes in population size
and stability were computed for each of the climate models as
(future − historical) /historical × 100% and plotted without outliers in
Fig. 4. Statistically significant changes in population abundance and
stability between the historical and future periods were identified via
the Mann–Whitney U-test, with the eight models as replicates.

Extinction probability was quantified as the proportion of ensem-
ble simulations for which the population declined to zero during a
50 yr simulation. Changes in extinction probability were calculated as
the difference between future and historical extinction probabilities.
Statistically significant changes in extinction probability were identi-
fied on a regional basis via the Mann–Whitney U-test.

Reporting summary
Further information on research design is available in the Nature
Research Reporting Summary linked to this article.

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 11

Data availability
The CMIP6 simulation data used in this paper are available via the data
portal https://esgf-node.llnl.gov/search/cmip6/. The ecology data are
available for download at https://doi.org/10.1073/pnas.0709472105.

Code availability
The code can be accessed on GitHub at https://github.com/KateDuffy/
climate-change-ecology54.

References
43. Gouhier, T. C., Grinsted, A. & Simko, V. R package biwavelet:
Conduct Univariate and Bivariate Wavelet Analyses (Version
0.20.17).

44. Torrence, C. & Compo, G. P. A practical guide to wavelet analysis.

Bull. Am. Meteorol. Soc. 79, 61–78 (1998).

45. Grinsted, A., Moore, J. C. & Jevrejeva, S. Application of the cross

wavelet transform and wavelet coherence to geophysical time
series. Nonlinear Process. Geophys. 11, 561–566 (2004).
46. Liu, Y., Liang, X. S. & Weisberg, R. H. Rectification of the bias

in the wavelet power spectrum. J. Atmos. Ocean. Technol. 24,
2093–2102 (2007).

47. Cazelles, B. et al. Wavelet analysis of ecological time series.

Oecologia 156, 287–304 (2008).

48. Maraun, D. & Kurths, J. Cross wavelet analysis: significance testing
and pitfalls. Nonlinear Process. Geophys. 11, 505–514 (2004).
49. Huey, R. B. & Stevenson, R. D. Integrating thermal physiology and
ecology of ectotherms: a discussion of approaches. Am. Zool. 19,
357–366 (1979).

50. Mallet, J. The struggle for existence: how the notion of carrying
capacity, K, obscures the links between demography, Darwinian
evolution, and speciation. Evol. Ecol. Res. 14, 627–665 (2012).
51. Vasseur, D. A. Theoretical Ecology: Concepts and Applications
(eds. McCann, K. S. & Gellner, G.) 243–262 (Oxford Univ. Press,
2020).

52. Estay, S. A., Clavijo-Baquet, S., Lima, M. & Bozinovic, F. Beyond
average: an experimental test of temperature variability on the
population dynamics of Tribolium confusum. Popul. Ecol. 53,
53–58 (2011).

53. Bozinovic, F. et al. The mean and variance of environmental

54. Duffy, K. climate-change-ecology (GitHub, 2022); https://github.

com/KateDuffy/climate-change-ecology

Acknowledgements
This work was primarily supported by the National Science
Foundation (NSF) grant CCF-1442728 while K.D. was a PhD
student at the SDS Lab in Northeastern University. Additional
support was provided for K.D. and A.R.G. by NSF SES-1735505
and for T.C.G. by NSF OCE-2048894. We acknowledge the
background support from a previous NSF Expeditions in
Computing grant (award no. 1029711) and an ongoing DOD
Strategic Environmental Research and Development Program
funding (no. RC20-1183). K.D. and A.R.G. acknowledge support
from the NASA Ames Research Center.

Author contributions
K.D., T.C.G. and A.R.G. conceived, designed and refined the project.
K.D. performed the data analysis and modelling. K.D., T.C.G. and A.R.G.
interpreted the results. K.D. wrote the manuscript with contributions
from T.C.G. and A.R.G.

Competing interests
The authors declare no competing interests.

Additional information
Extended data is available for this paper at
https://doi.org/10.1038/s41558-022-01490-7.

Supplementary information The online version
contains supplementary material available at
https://doi.org/10.1038/s41558-022-01490-7.

Correspondence and requests for materials should be addressed to
Kate Duffy.

Peer review information Nature Climate Change thanks David Vasseur
and the other, anonymous, reviewer(s) for their contribution to the
peer review of this work.

temperature interact to determine physiological tolerance and
fitness. Physiol. Biochem. Zool. 84, 543–552 (2011).

Reprints and permissions information is available at
www.nature.com/reprints.

Nature Climate Change

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 12

Extended Data Fig. 1 | Temperature variation at multiple timescales
contributes to trends in overall variance. a-b Temporal trends in the power
of variation at sub-annual to annual periodicities (3-days to 2-years) (a) and
multi-annual periodicities (2-30 years) (b). Trends represent the slope obtained

by regressing wavelet power at each geographical location against time.
Countervailing trends are found in the Arctic, where the power of short term,
high-frequency fluctuations is decreasing and the power of more persistent,
low-frequency fluctuations is increasing.

Nature Climate Change

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 13

Extended Data Fig. 2 | Drivers of changes in stability (analysis includes
both pre- and post-extinction period). Kernel density plots illustrate the
relationships between population mean and population standard deviation
in the historical period and the future climate change period. The grey 1:1 line
divides the more stable regime (high-mean/low-variance; below line), and the

less stable regime (low-mean/high-variance; above line). Bimodal distributions
emerge in the extra-tropics, with some species at low abundance and standard
deviation, and a larger cluster of species at high abundance and standard
deviation.

Nature Climate Change

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 14

Extended Data Fig. 3 | Drivers of changes in stability (analysis only includes pre-extinction period). When only pre-extinction dynamics are analyzed, significant
changes in population abundance persist in the extra-tropics; changes in population standard deviation become significant for NHEX and TROP and remain non-
significant for SHEX.

Nature Climate Change

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 15

Extended Data Fig. 4 | Scaling of the intrinsic growth rate has moderate
effects on temperature-driven impacts on population stability and
extinction risk. Results exhibited limited sensitivity to the choice of smaller
(scaling factor = 0.1; a,b) and larger (scaling factor = 10.0; e,f) intrinsic growth

rates. Although larger growth rates were more strongly associated with deceased
stability and increased extinction risk than smaller growth rates, the latitudinal
patterns and effect sizes were consistent with the changes in population stability,
c, and extinction probability, d, observed under normalized growth rates.

Nature Climate Change

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 16

Extended Data Fig. 5 | Temperature-driven effects on population stability and extinction risk are robust to the degree of population self-regulation. Results
exhibited limited sensitivity to strong (α = 1; Fig. 4) and weak (α = 0.1; above) self-regulation in the form of crowding effects. Latitudinal patterns and effect sizes were
consistent for changes in population stability, a, and extinction probability, b.

Nature Climate Change

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 17

Extended Data Fig. 6 | Patterns of species risk are sensitive to treatment of
high-temperature performance. Change in the average fitness (growth rate)
driven by daily temperature between the historical period (1950-2000) and
future period (2050-2100) when allowing negative growth rates above CTmax
as in Vasseur et al. 2014 (a) and when performance values are bounded by 0

above CTmax as in Deutsch et al. 2008 (b). Change in the average fitness (growth
rate) driven by monthly mean temperature between the historical period and
future period when allowing negative growth rates above CTmax (c) and when
performance values are bounded by 0 above CTmax (d).

Nature Climate Change

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 18

Extended Data Fig. 7 | Increased stability is negatively related to extinction
probability. Regression relationships in our simulations are presented a, when
considering only the pre-extinction time period and b, when taking into account

the full 50-year periods. Regardless of largely positive (b) or mixed (a) changes in
stability, there is generally a weak but significant negative relationship between
stability and extinction probability globally (p-value < 0.05).

Nature Climate Change

Articlehttps://doi.org/10.1038/s41558-022-01490-7

Page 19

Corresponding author(s):

Kate Duffy

Last updated by author(s): Aug 23, 2022

Reporting Summary

Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency
in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed

The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement

A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly

The statistical test(s) used AND whether they are one- or two-sided
Only common tests should be described solely by name; describe more complex techniques in the Methods section.

A description of all covariates tested

A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons

A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient)
AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals)

For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted
Give P values as exact values whenever suitable.

For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings

For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes

Estimates of effect sizes (e.g. Cohen's d, Pearson's r), indicating how they were calculated

Our web collection on statistics for biologists contains articles on many of the points above.

Software and code

Policy information about availability of computer code

Data collection

No software was used for data collection in this study.

Data analysis

Open source codes were used to analyze the data, including R package quantreg v5.88 and R package biwavelet v0.02.21. Custom codes for
analysis are available in our public GitHub repository at https://github.com/KateDuffy/climate-change-ecology

For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and
reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information.

Data

Policy information about availability of data
All manuscripts must include a data availability statement. This statement should provide the following information, where applicable:

- Accession codes, unique identifiers, or web links for publicly available datasets
- A description of any restrictions on data availability
- For clinical datasets or third party data, please ensure that the statement adheres to our policy

The CMIP6 simulation data used in this paper is available via the data portal https://esgf-node.llnl.gov/search/cmip6/. The ecology data is available for download at
https://doi.org/10.1073/pnas.0709472105.

n
a
t
u
r
e
p
o
r
t
f
o

l
i

r
e
p
o
r
t
i
n
g
s
u
m
m
a
r
y

M
a
r
c
h
2
0
2
1

Page 20

Human research participants

Policy information about studies involving human research participants and Sex and Gender in Research.

Reporting on sex and gender

Not applicable

Population characteristics

Not applicable

Recruitment

Ethics oversight

Not applicable

Note that full information on the approval of the study protocol must also be provided in the manuscript.

Field-specific reporting

Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.

Life sciences

Behavioural & social sciences

Ecological, evolutionary & environmental sciences

For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf

Ecological, evolutionary & environmental sciences study design

All studies must disclose on these points even when the disclosure is negative.

Study description

Research sample

Using simulations from the latest generation of earth system models, we applied quantile, spectral, and wavelet analyses of
temperature projections to study the way that temperatures will increase over time. We integrated temperature projections into
empirically-parameterized mathematical models that simulate the dynamical and cumulative effects of thermal stress on the
performance of 38 global ectotherm species.

We obtained experimentally derived thermal tolerance parameters for a set of terrestrial ectotherms of the phylum arthropoda
(n = 38) published by Deutsch et al. (2008). Deutsch gathered data from 31 previously published thermal performance studies, which
were published between 1974 and 2003 based on a collection of insects from 35 different locations. For each species, experimental
intrinsic growth rates at multiple temperatures were used to fit a thermal performance curve yielding least-squares estimates of the
the curve's parameters. This dataset was selected to represent a globally distributed set of arthropod species.

Sampling strategy

Not applicable to study. We analyzed existing data as described in the data availability statement.

Data collection

Not applicable to study. We analyzed existing data as described in the data availability statement.

Timing and spatial scale Not applicable to study. We analyzed existing data as described in the data availability statement.

Data exclusions

No data were excluded from the analysis.

Reproducibility

Our study consists of fully reproducible data analysis.

Randomization

Not applicable.

Blinding

Not relevant to study.

Did the study involve field work?

Yes

Reporting for specific materials, systems and methods

We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material,
system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.

n
a
t
u
r
e
p
o
r
t
f
o

l
i

r
e
p
o
r
t
i
n
g
s
u
m
m
a
r
y

M
a
r
c
h
2
0
2
1

Page 21

Materials & experimental systems
n/a Involved in the study

Methods
n/a Involved in the study

Antibodies

Eukaryotic cell lines

Palaeontology and archaeology

Animals and other organisms

Clinical data

Dual use research of concern

ChIP-seq

Flow cytometry

MRI-based neuroimaging

n
a
t
u
r
e
p
o
r
t
f
o

l
i

r
e
p
o
r
t
i
n
g
s
u
m
m
a
r
y

M
a
r
c
h
2
0
2
1

Page: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21